******SET YOUR DIRECTORY to ...\replication



clear all


/////////////////////////////////////////////////////////////////////////////////////////////
///Proportion with water/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////

***Map only intersection radius 20 and places
shp2dta using intermediatedata\inter_radius20_places_naec, data("map-­-attr_pl.dta") coord("map-­-coord_pl.dta") ///
 genid(stid)  gencentroids(cc) replace
 
use intermediatedata\map-­-attr_pl.dta, clear

gen fipsplace_00=string(fipsplace) if fipsplace>=1000000
replace fipsplace_00="0"+string(fipsplace) if fipsplace<1000000

*Negative land value for Anchorage
replace ALAND00=abs(ALAND00) if fipsplace_00=="0203000"

gen area_water_places_sqkm=(AREA*(AWATER00/(AWATER00+ALAND00)))/1000000

collapse (sum) area_water_places_sqkm, by(fipsplace_00)

save rawdata\inwatera_usa\data_undevelopable_m.dta, replace

***Map only intersection radius 20 and coastal water (without intersection with places and inland water)
shp2dta using inter_radius20_coastwater_naec_nopliw, data("map-­-attr_cw.dta") coord("map-­-coord_cw.dta") ///
 genid(stid)  gencentroids(cc) replace
 
use intermediatedata\map-­-attr_cw.dta, clear

gen fipsplace_00=string(fipsplace) if fipsplace>=1000000
replace fipsplace_00="0"+string(fipsplace) if fipsplace<1000000

gen area_water_cw_sqkm=(AREA/1000000)

collapse (sum) area_water_cw_sqkm, by(fipsplace_00)

merge 1:1 fipsplace_00 using rawdata\inwatera_usa\data_undevelopable_m.dta
drop _merge 

save rawdata\inwatera_usa\data_undevelopable_m.dta, replace

***Map only intersection radius 20 and inland water (without intersection with places)
shp2dta using inter_radius20_inwater_naec_nopl, data("map-­-attr_iw.dta") coord("map-­-coord_iw.dta") ///
 genid(stid)  gencentroids(cc) replace
 
use intermediatedata\map-­-attr_iw.dta, clear

gen fipsplace_00=string(fipsplace) if fipsplace>=1000000
replace fipsplace_00="0"+string(fipsplace) if fipsplace<1000000

gen area_water_iw_sqkm=(AREA/1000000)

collapse (sum) area_water_iw_sqkm, by(fipsplace_00)

merge 1:1 fipsplace_00 using rawdata\inwatera_usa\data_undevelopable_m.dta
drop _merge 

***Generate total area water 
replace area_water_iw_sqkm=0 if area_water_iw_sqkm==.
replace area_water_cw_sqkm=0 if area_water_cw_sqkm==.
replace area_water_places_sqkm=0 if area_water_places_sqkm==.
gen area_water_sqkm=area_water_iw_sqkm+ area_water_cw_sqkm+ area_water_places_sqkm
gen prop_area_water=area_water_sqkm/(20*20*3.14)

save rawdata\inwatera_usa\data_undevelopable_m.dta, replace

/////////////////////////////////////////////////////////////////////////////////////////////
///Proportion with slope/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////
///Convert  shapefile  to  Stata  attribute  and  coordinate  datasets
shp2dta using inter_radius20_statsgo_naec_nocwiw, data("map-­-attr_slope.dta") coord("map-­-coord_slope.dta") ///
 genid(stid)  gencentroids(cc) replace

use intermediatedata\map-­-attr_slope.dta, clear

destring MUKEY, gen(mukey) force
drop MUKEY
gen fipsplace_00=string(fipsplace) if fipsplace>=1000000
replace fipsplace_00="0"+string(fipsplace) if fipsplace<1000000
keep mukey fipsplace_00 AREA

save intermediatedata\map-­-attr_slope_m.dta, replace

///Merge with mapunit.txt
insheet using "rawdata\wss_gsmsoil_US_[2006-07-06]\tabular\mapunit.txt", clear delimiter("|")
rename v1 musym
rename v2 muname
rename v3 nukind
rename v4 mustatus
rename v5 muacres
rename v6 mapunitflfw_l
rename v7 mapunitflfw_r
rename v8 mapunitflfw_h
rename v9 mapunitpfa_l
rename v10 mapunitpfa_r
rename v11 mapunitpfa_h
rename v12 farmlndcl
rename v13 muhelcl
rename v14 muwathelcl
rename v15 muwndhelcl
rename v16 interpfocus
rename v17 invesintens
rename v18 iacornsr
rename v19 nhiforsoigrp
rename v20 nhspiagr
rename v21 vtsepticsyscl
rename v22 mucertstat
rename v23 lkey
rename v24 mukey
label var mukey "Map unit ID"
label var muname "Map unit name"

keep mukey muname

merge 1:m mukey using intermediatedata\map-­-attr_slope_m.dta
drop if _merge==1
drop _merge

save intermediatedata\map-­-attr_slope_m.dta, replace

///Merge comp.txt with chorizon.txt
insheet using "rawdata\wss_gsmsoil_US_[2006-07-06]\tabular\comp.txt", clear delimiter("|")
rename v2 comppct_r
rename v4 compname
rename v10 slope_r
rename v26 elev_r
rename v39 map_r
rename v108 mukey
rename v109 cokey
label var comppct_r "Percentage of the component of the mapunit"
label var compname "Component name"
label var slope_r "Slope RV"
label var elev_r "Elevation RV"
label var map_r "Arithmetic average of the total annual (liquid) precipitation taken over 1961-1990"
label var mukey "Map unit ID"
label var cokey "Component ID"

drop v*

///Collapse at map unit level
egen slope_wtm_r=wtmean(slope_r),by(mukey) weight(comppct_r)

collapse (mean) slope_wtm_r,by(mukey)

///Merge with map units
merge 1:m mukey using intermediatedata\map-­-attr_slope_m.dta
drop if _merge==1
drop _merge

***Generate total area slope
drop if slope_wtm_r<15
gen area_slope_sqkm=(AREA/1000000)

collapse (sum) area_slope_sqkm, by(fipsplace_00)

gen prop_area_slope=area_slope_sqkm/(20*20*3.14)

merge 1:1 fipsplace_00 using rawdata\inwatera_usa\data_undevelopable_m.dta
drop _merge 

replace prop_area_slope=0 if prop_area_slope==.

gen prop_undevelopable=prop_area_slope+prop_area_water

destring fipsplace_00, force replace

save rawdata\inwatera_usa\data_undevelopable.dta, replace
